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SUMMARY 


An effort to develop a three-dimensional ice accretion modeling method has been initiated. This 
first step towards creation of a complete aircraft icing simulation code builds on previously developed 
methods for calculating three-dimensional flowfields and particle trajectories combined with a two- 
dimensional ice accretion calculation along coordinate locations corresponding to streamlines. This 
work is intended as a demonstration of the types of calculations necessary to predict a three- 
dimensional ice accretion. Results of calculations using the 3D method for a MS-317 swept wing 
geometry are projected onto a 2D plane normal to the wing leading edge and compared to 2D results 
for the same geometry. It is anticipated that many modifications will be made to this approach, 
however this effort will lay the groundwork for future modeling efforts. Results indicate that the 
flowfield over the surface and the particle trajectories differed for the two calculations. This led to 
lower collection efficiencies, convective heat transfer coefficients, freezing fractions, and ultimately ice 
accumulation for the 3D calculation. 
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I. Introduction 

Ice accretion modeling has received considerable attention during the past few years. The 
National Aircraft Icing Technology Plan 1 , drafted in 1986 under the direction of several federal 
agencies, calls for the creation of a complete aircraft icing simulation during the coming decade. In 
support of this plan several activities related to the creation of a 2D ice accretion modeling and 
performance evaluation capability have been undertaken at the NASA Lewis Research Center. A — D 
ice accretion code, LEWICE, was developed in 1983 by the University of Dayton Research Institute 2 
and later modified by Ruff 3 . This code is divided into three major parts; a potential flow calculation, 
a particle trajectory calculation, and an ice growth calculation. The ability of LEWICE to predict 
airfoil ice accretions is quite good for rime ice growths. Glaze ice predictions, although acceptable, 
have exhibited some deficiencies which have been documented by several authors . The calculation 
of performance degradation due to icing is an equally important aspect of the icing analysis problem. 
Over the past several years 2D calculations for iced airfoils have been performed using the interactive 
boundary layer method 2 and using Navier-Stokes codes . Recently, Cebeci has combined his 
interactive boundary layer method with LEWICE. This has resulted in a code which has the potential 

to accurately determine the flowfield and particle trajectories while also providing the performance 
degradation information during a single calculation. 
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With these activities firmly underway, it is appropriate to consider the methodology for 
development of a three-dimensional ice accretion analysis capability. Some of the issues associated 
with a 3D icing calculation are the same as those for 2D. On the other hand, some issues only become 
evident when considering a full 3D calculation. These may include; variations in ice growth on both a 
large and small scale across a wing span, runback in the spanwise direction, non-uniformities in the 
impinging cloud along the span, and particle trajectories in complex geometric regions. 

The current effort is directed at taking existing computer codes and developing a strategy for 
performing a 3D ice accretion and performance degradation calculation for a swept wing configuration. 
This activity will serve as a demonstration calculation and as a focus for future development of the 
approach. The computer codes used in this effort consist of a 3D panel code, a particle trajectory code, 
and the ice accretion subroutines of the LEWICE code. The performance calculations will use a 3D 
Navier-Stokes code, ARC3D. Presently, the calculations are performed separately. In the future, 
these calculations will be performed in a single code thus requiring less user interaction than is 
currently necessary. This report describes the ice accretion calculation as applied to an infinite swept 
wing configuration. 

The wing geometry chosen for this calculation is that of an MS-317 airfoil with a 30° sweep angle, 
as shown in Figure I. The MS-317 airfoil model represents a typical medium-speed wing section. This 
geometry was chosen due to the availability of experimental information from ice accretion tests 
currently being carried out in the NASA Lewis Icing Research Tunnel. Comparisons of the predicted 
ice accretion profiles to those of the experiment will be completed as a follow-on to this work. 

II. Code Description 

A 3D ice accretion calculation for a complete aircraft is a major step forward from current 
capabilities. The present activity is directed at extending current 2D capabilities to simple 3D 
geometries by assuming constant spanwise gradients. For this work, an infinite swept wing of constant 
geometry in the spanwise direction is used. Such a geometry produces similar streamlines along the 
span. As a result of this streamline similarity, it is assumed that there is no net transport of water 
across the control volume in the spanwise direction. The control volume analysis of the ice accretion 
can thus be considered two dimensional and the 2D LEWICE routines can be used. This also allows a 
single evaluation of the ice growth for any streamline to represent the ice accretion behavior for any 
location on the wing. Since the resulting ice shape geometry will not vary along the span, re-paneling 
of the iced wing for subsequent ice accretion calculations is simplified. 

An infinite swept wing also avoids consideration of ice growth near the tip of the wing or at the 
wing-body junction. At these locations, the flowfield is quite irregular and major changes to the 
approach currently being employed would be required. It is anticipated that future work in 3D ice 
accretion prediction would be directed at addressing these regions of the flowfield. Currently, the tools 
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necessary for evaluation of these complex flowfields are being developed by Sankar 10 and Caruso 11 . 
Additional work will be required to create a 3D control volume analysis for the ice accretion 
calculation. As these codes become available, they will be incorporated into the 3D ice accretion 
analysis methodology. 

II- 1. 3D Panel Code and Particle Trajectory Calculation 

The flowfield and particle trajectory calculations are extensions of the 2D analysis methods used 
in LEWICE. The three-dimensional flowfield calculation uses a Hess-Smith panel code 12 . The code 
can accomodate lifting and non-lifting geometries or combinations, thus allowing simulation of complex 
surfaces such as a complete aircraft (Figure 2). A Prandtl-Glauert correction is used to allow 
compressible flow calculations. The code can also handle leaking panels to simulate inlets or 
instrument orifices. The code provides velocity information at any point in the flowfield away from the 
surface. This type of code uses very little computer time compared to a viscous flow code and yet 
provides sufficient flowfield information for subsequent trajectory and energy balance calculations under 
most circumstances. 

The trajectory calculation is based on the computer code developed by Norment 13 with one 
additional feature. The code solves a force balance equation on a single particle and then determines 
the new particle position using an Adams-type predictor-corrector algorithm developed by Krogh 14 . 
The added feature is the ability to calculate the local collection efficiency from the impacting particles. 
The code generates an array of impingement points for each region of interest. This is done by 
releasing an array of particles which impact the region of interest (Figure 3). The collection efficiency 
is the ratio of the particle flux at the target point to the free-stream particle flux. If a known group of 
particles is tracked from their release location to the target location, then the collection efficiency can 
be determined by the ratio of the target area to the release area. The collection efficiency can then be 
determined from the following relationship, 

/?(xc(ij)) = A 0 (ij)/ A m (i j) (1) 

Once the collection efficiency is determined for the surface, an array of surface locations and 
corresponding fi values is created for streamlines on the surface. This information will be used in the 
ice accretion calculation to determine ice growth along the streamline paths. The streamlines are 
calculated using a 4th order Runge-Kutta integration scheme. The streamline is carried forward from 
the stagnation region for both the upper and lower surface at the region of interest. 

A linear interpolation scheme is used to determine the collection efficiency along the streamline 
from the matrix of (3 values generated in the trajectory calculation. The array of (3 values is searched 
to find the surface cell in which the streamline point resides (Figure 4). The (3 value at the streamline 
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point is then determined from a weighted average of the local collection efficiencies at the four 
surrounding cell points. The interpolation expression is, 


P{x 3 ) ~ /?(x c (i+ld + l))-u-v + /?(x c (ij + l))-u-(l-v) + /?(x c (ij))-(l-u)-(l-v) 
+ /?(x c (i+lj))-(l-u)-v 


where, 


and 


u = X, • Xj 
v = X, • Xj 


( 2 ) 

( 3 ) 

( 4 ) 


represent weighting factors in the interpolation based on the geometry of the cell. 


IT-2. Ice Accretion Calculation 

The ice accretion calculation is performed as a 2D strip analysis along streamlines calculated by 
the 3D panel code. The results of the previous calculation are written into a file containing the x,y,z 
coordinates of the streamline, the surface distance between points, the x,y,z components of the unit 
normal vector for the panel at the streamline coordinate locations, the velocities, and the local 
collection efficiencies, f3. The subroutines from the LEWICE code which determine the boundary layer 
values, the control volume energy balance and the resulting ice growth are used to determine the ice 
shape at a given spanwise location. A brief outline of the LEWICE calculation follows. Further details 
can be found in reference 3. 

The ice accretion calculation consists of a control volume mass and energy balance. The 2D strip 
analysis allows use of the assumption that there are no spanwise fluxes into the control volume, which 
is consistant with the infinite swept wing flowfield assumption. A depiction of the control volume for 
the process is shown in Figure 5. The As value is based on distance along the streamline. 

The mass balance equation is, 


m c + rii r . n = m e + m r<>u( + ihj (5) 

The impinging mass flux, rii c , is determined from the (5 distribution. The mass flux into the control 
volume, m r . , is any water remaining in the liquid state after evaluation of an upstream control 
volume. This mass flux is called runback. The mass flux out of the control volume due to 
evaporation, m e , and to freezing, im, are determined from the energy balance calculations. The mass 
flux out of the control volume, m r , is used as input for the downstream control volume. The 
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freezing fraction is a term used to describe the fraction of incoming liquid that freezes within the 
control volume. It is given by the expression, 


/ = 


m : 


m c + m r»n 


( 6 ) 


The energy balance equation is, 


qc + Qr #n — q c on d + Qevap + <\conv + qr ou< + qj (7) 

The values of q c and q r are obtained from mass flux values and the internal energy levels of 
the incoming water. Integral boundary layer expressions are used to determine the heat flux due to 
convection. The inviscid flowfield values are used as boundary values for the viscous calculation. The 
boundary layer calculation determines the location of transition from laminar to turbulent flow and 
then uses an appropriate value of the convection heat transfer coefficient. The heat transfer value is 
then obtained by use of Newton’s law of cooling, 

<‘|ccnv = ll A (T, - Too) ( 8 ) 

The evaporative heat flux is also obtained from consideration of the inviscid flowfield conditions using 
the expression, 


q tvap 


— (Cp 


T, -f Lt^As 


( 9 ) 


where the evaporative mass transfer rate, ih e , is determined from local pressure and temperature 
conditions. Further details on this calculation are found in reference 3. The heat flux due to 
conduction is a user specified value. This leaves two equations for the two unknowns, rh r< ^ and riij. 

The mass of ice that freezes is used to determine the thickness of the resulting ice layer. Since the 
analysis is locally 2D, the ice thickness can be found from the relation, 


^tce 


m. At 
1 


Pic 


As W 


( 10 ) 


The thickness is then considered to be uniform over tire entire panel for subsequent flowfield 
calculations. The new coordinates for the panel are obtained from the relation, 


x, = x, + d lce x< 


( 11 ) 


where x, is the *th coordinate of the center of the panel and x^ is the ith component of the unit normal 
vector for the panel. 


6 


As the ice thickness increases, there is the possibility that ice segments will intersect and thus this 
must be accounted for in the determination of the new geometry. Since this is a strip analysis, the 
spanwise thickness does not vary. Therefore, the possibility of ice growth intersection is limited to the 
normal and chordwise directions. In that case, the line segments corresponding to the top of every 
other panel are examined for intersection. If intersection is determined to occur, then a new panel is 
formed with its center halfway between the two old panels. This requires determination of the 
coordinates of the new panel and renumbering of the panels. This information is then used in 
subsequent potential flow calculations. 


III. Results and Discussion 

Prior to the development of modern computational techniques, Dorsch and Brun 15 devised a 
method for determining cloud-droplet impingement on swept wings. Their approach is based on the 
argument that droplet trajectories for an infinite swept wing can be calculated in the same manner as 
for a 2D plane normal to the leading edge of the wing, using the component of the free-stream velocity 
in that plane. Thus, the impingement limits, local collection efficiencies, and total water collected are 
calculated for the normal plane. The results are then considered to be equivalent at all normal planes 
along the span. 

An examination of that approach can be undertaken by comparing results of the current method 
to a 2D LEWICE calculation for a normal plane of the swept wing configuration. In order to compare 
the swept to unswept condition, the approach velocity of the 2D calculation was taken as the 
component of the swept wing velocity in the normal plane. If 7 is the sweep angle, then the 
appropriate 2D approach velocity is found using the relation, 

V OOti — V OO s cos 7 ( 12 ) 


In order to understand the forces exerted on the incoming particles, it is illustrative to compare 
the relative strengths of the pressure fields for the two cases to the inertia of the incoming water 
droplets. The inertia of the droplets is characterized by the expression, 

'iPwVoo 2 (13) 


An alternate pressure coefficient based on the droplet inertia would have the form, 


Cn = 


p — Poo 


Jp y w Voo 2 


(14) 


7 


This term is simply related to the standard pressure coefficient through the ratio of water to air 
densities. Thus, the standard Cp values based on the density of air and the approach velocity can be 
used to evaluate the forces exerted on the water droplets. These are plotted in Figure 6 for the two 
cases considered. 

The 3D potential flow calculation for the MS-317 wing section with a 30° sweep angle yielded 
pressure coefficient profiles in a plane normal to the leading edge as shown in Figure 6. The Cp 
distribution is contrasted to that of a 2D airfoil section in the figure to show the influence of the sweep 
angle. The swept wing Cp distribution has somewhat higher pressure values on the suction surface. 
Based on the discussion above, this indicates that the swept wing suction pressure has less influence on 
the incoming particle trajectories than the 2D suction pressure. This difference causes a shifting of the 
local collection efficiency distribution, as discussed below. The differences in flowfield characteristics 
also have an impact on the convective heat transfer in the region of ice accretion. 

Three-dimensional particle trajectory calculations from this flowfield are shown in Figure 7. The 
trajectories indicate particle motion in the spanwise direction in the region just in front of the wing. 
Apparently the influence of the spanwise velocities is not felt until the particles are close to the surface. 
The free-stream velocity, sweep angle, and particle size all play a role in the extent of spanwise travel. 
The particles thus impact the surface at some angle with respect to the chordwise direction. The effect 
of this 3D particle trajectory is evident in the local collection efficiency plots, as shown in Figure 8. 
Additionally, more complex geometries would affect particle trajectories to a greater degree than the 
simple geometry examined here. 

Along the streamline, the local collection efficiencies are determined according to Eq. 2. For 
comparison purposes, a plot of (3 vs. distance from the stagnation point in a plane normal to the 
leading edge of the wing, along with 2D results for the same geometry, are shown in Figure 8. The 3D 
ft curve has impingement limits at a distance somewhat closer to the stagnation point than the 2D 
case. Additionally, the 3D results indicate a somewhat lower (3 near the stagnation region. This 
behavior is most noticeable for the suction surface (i.e. positive s values). The higher pressure levels on 
the suction surface of the 3D calculation result in a relative shift of the particle trajectories to the 
pressure surface compared to the 2D case. This combined with spanwise motion of the particles 
described above account for the different (3 distributions. 

The convective heat flux term has a major influence on the control volume energy balance. The 
velocity distribution over the surface and the resulting convective heat transfer coefficient are shown in 
Figures 9 and 10 respectively. Since the spanwise velocity does play a role in the convective heat 
transfer of the energy balance, the total free-stream velocity was used in this portion of the 2D 
calculation. The velocity distributions are somewhat different for the two calculations. The higher 
pressure of the 3D calculation is reflected in the lower velocity levels. This results in lower heat 
transfer coefficient values for that surface. The greatest differences are in the stagnation region. The 
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milder pressure-side inviscid velocity gradients from the 3D calculation result in a lower heat transfer 
coefficient value. This suggests the possibility of more runback water on the pressure side for the 3D 
calculation. Additionally, the 2D calculation produces smoother distributions. Examination of 
interpolation routines and surface panel distribution is required. 

These changes in the heat transfer distribution on the surface combined with the altered j3 
distribution lead to differences in the surface temperature and freezing fraction, as shown in Figures 11 
and 12 respectively. Both figures show the shifting of the region of ice accumulation to the pressure 
side of the wing. The surface temperature plots show a small region of transition from an ice/water 
interface (indicated by the regions where T is constant at 0°C) to a completely dry, iced surface 
(indicated by the regions where T is constant at approximately -16°C) for the 2D calculation. The 3D 
calculation indicates a dry iced surface over the entire accretion area. The large gradient in surface 
temperature at s = 0.024 for the 3D calculation may be due to the sharp spikes in the heat transfer 
coefficient distribution. The freezing fraction plots indicate a small region of wet ice growth for the 2D 
calculation and a completely dry ice growth for the 3D calculation. The heat transfer coefficient in this 
region does not drop as low for the 3D calculation as it does in the 2D calculation. This results in 
sufficient energy loss to provide for freezing of all the incoming water. 

The ice accretion results for this wing are shown in Figure 13. Comparison to the 2D results 
indicate the influence of the altered collection efficiency distribution and of the 3D flowfield. The ice 
shape obtained with the 3D calculation has less mass and extends over a smaller portion of the surface. 
In this calculation, that shape would then be considered to be uniform along the span. 

IV. Conclusion 

An approach to calculate the ice accretion for an infinitely swept wing has been proposed. This 
method employs existing technology to couple the evaluation of 3D particle trajectories with a 2D strip 
analysis to predict ice buildup along surface streamlines. This reduces the computational resources that 
would be required for predicting the full 3D flowfield and uses only a small increase in resources over 
the current 2D LEWICE code. This approach could be easily extended to other simple 3D geometries 
such as axisymmetric engine inlets or rotor blades. 

Comparisons to a 2D analysis were made for a MS-317 infinite swept wing geometry. Results 
indicate that the flowfield over the surface and the particle trajectories differed for the two calculations. 
This led to lower collection efficiencies, convective heat transfer coefficients, freezing fractions, and 
ultimately ice accumulation for the 3D calculation. These results also indicate that improvements can 
be made to the previously used method for evaluation of swept wing cloud-droplet impingement. The 
present calculations include the influence of sweep on pressure levels relative to droplet inertia. 

Additional work is required to determine the number of panels necessary to accurately model the 
flowfield and resulting collection efficiency distribution for the swept wing. Studies of several swept 
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wing geometries are necessary to confirm the current results. Comparison to actual ice accretions for a 
swept wing geometry are also required to provide confidence in the calculation procedure. Future 
studies of the effects of finite span and wing/body junctions will be incorporated to give a prediction of 
the ice accretion distribution on complete wing geometries. 
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Figure 1 - MS-317 Swept Wing Profile 
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Finely resolved version 



Coarsely resolved version 


Digital description of the DeHaviiland Twin Otter. 
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Figure 5. Control Volume for Ice Accretion Calculation 
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PRESSURE COEFFICIENT DISTRIBUTION 





Figure 6 - Pressure Coefficient Distribution 
for MS-317 Profile: 2D and 3D Results 



Figure 7 - Particle Trajectories and Streamline 
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Figure 8. Local Collection Efficiency Distribution in a 
Plane Normal to the Leading Edge for a MS-317 
Swept Wing : 2D and 3D Calculations 
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MS-317 : 2D CALCULATION 
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Figure 9. Local Inviscid Velocity Distribution in a Plane 
Normal to the Leading Edge for a MS-317 Swept 
Wing : 2D and 3D Calculations 
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Figure 10. Convection Heat Transfer Coefficients in a 

Plane Normal to the Leading Edge for a MS-317 
Swept Wing : 2D and 3D Calculations 
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Figure 11. Surface Temperature Near Stagnation Region of 
MS-317 Swept Wing : 2D and 3D Calculations 
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Figure 12. Freezing Fraction Near Stagnation Region of 
MS-317 Swept Wing ; 2D and 3D Calculations 
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Figure 13. Ice Shape Profiles Near Leading Edge of 

MS-317 Swept Wing ; 2D and 3D Calculations 


24 





!x§^csan d Report Documentation Page 

Space Administration 

j 1. Report No. TM^ / 2 Government Accession No. 

i AIAA-90-0756 

3. Recipient’s Catalog No. 

4. Title and Subtitle 

Swept Wing Ice Accretion Modeling 

5. Report Date 

6. Performing Organization Code 

7 Author(s) 

Mark G. Fotapczuk and Colin S. Bidwell 

8. Performing Organization Report No 

E-5238 

10. Work Unit No. 

505-68-11 

9. Performing Organization Name and Address 

National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio 44135-3191 

11. Contract or Grant No. 

13. Type of Report and Period Covered 
Technical Memorandum 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, D.C. 20546-0001 

14. Sponsoring Agency Code 

15. Supplementary Notes 

Prepared for the 28th Aerospace Sciences Meeting sponsored by the American Institute of Aeronautics and 
Astronautics, Reno, Nevada, January 8-11, 1990. 


16. Abstract 


An effort to develop a three-dimensional ice accretion modeling method has been initiated. This first step towards 
creation of a complete aircraft icing simulation code builds on previously developed methods for calculating 
three-dimensional flowfields and particle trajectories combined with a two-dimensional ice accretion calculation 
along coordinate locations corresponding to streamlines. This work is intended as a demonstration of the types of 
calculations necessary to predict a three-dimensional ice accretion. Results of calculations using the 3D method 
for a MS-317 swept wing geometry are projected onto a 2D plane normal to the wing leading edge and compared 
to 2D results for the same geometry. It is anticipated that many modifications will be made to this approach, 
however this effort will lay the groundwork for future modeling efforts. Results indicate that the flowfield over 
the surface and the particle trajectories differed for the two calculations. This led to lower collection efficiencies, 
convective heat transfer coefficients, freezing fractions, and ultimately ice accumulation for the 3D calculation. 


17. Key Words (Suggested by Author(s)) 
Aircraft icing 

Computational fluid dynamics 

18. Distribution Statement 

Unclassified — Unlimited 
Subject Category 02 

19. Security Classif. (of this report) 
Unclassified 

20. Security Classif. (of this page) 
Unclassified 

21 . No. of pages 

22. Price* 


NASA FORM 1626 OCT 86 


‘For sale by the National Technical Information Service, Springfield, Virginia 22161 





National Aeronautics and 
Space Administration 


SECOND CLASS MAIL 




Lewis Research Center address correction requested 

Cleveland. Ohio 44135 



Official Business 
Penalty lor Private Use $300 


Postage and Fees Paid 
National Aeronautics and 
Space Administration 
NASA -451 





ERRATA 


NASA Technical Memorandum 103114 


SWEPT WING ACCRETION MODELING 

Mark G. Potapczuk and Colin S. Bidwell 


Cover and Report Documentation Page: The NASA Technical 

Memorandum number should be changed from 102453 to 103114. 




